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We study the energy spectrum of a graphene bilayer in the presence of transverse electric and 
magnetic fields. We find that the resulting Landau levels exhibit a nonmonotonic dependence on 
the electric field, as well as numerous level crossings. This behavior is explained using quasiclas- 
sical quantization rules that properly take into account the pseudospin of the quasiparticles. The 
pseudospin generates the Berry phase, which leads to a shift in energy quantization and results in 
a pseudo-Zeeman effect. The latter depends on the electric field, alternates in sign among the two 
valleys, and also reduces the band gap. Analytic formulas for other pseudospin-related quantities, 
such as the anomalous Hall conductivity, are derived and compared with prior theoretical work. 
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I. INTRODUCTION 

The physics of monolayer and bilayer graphene has at- 
tracted much recent attention.^ A unique feature of bi- 
layer graphene (BLG) is its tunable band structure: the 
symmetric bilayer is gapless but when an interlayer po- 
tential difference U is induced, a band gap opens. The 
low-energy regions affected by the gap are situated at the 
Brillouin zone corners, e.g., points K — ±(47r/3a )5; 
henceforth referred to as K ± valleys, near which the band 
dispersion acquires a "Mexican hat" shape^(see Fig.[l|, 
where a — 2.46 A is the lattice constant for the under- 
lying triangular Bravais lattice. 

There have been interesting theoretical predictions 
that electron interactions can spontaneously generate 
layer polarization and a band gappEObut they have yet to 
be verified experimentally.^ The proven ways of creating 
an interlayer bias U include doping and gating. The lat- 
ter enables one to change U continuously, althoug h the 
dependence of U on the gate voltage is nontrivial! 11 * 12 ! 
In most of experimental studies of bilayer graphene a 
single gate electrode was usedP^HUl l n such devices the 
interlayer bias U and the induced electron density n vary 
concomitantly with the gate voltage. Sep arate control 
of U and n can be achieved with two gates I2212H Experi- 
ments with dual-gate devices^^MH] nave been reported 
recently. 

Another intriguing property of graphene is that its low- 
energy quasiparticles are endowed with a pseudospin- 1 
degree of freedom, associated with the sublattice struc- 
ture of each monolayer, whose dynamics is linked to 
their orbital motion. 1 When a quasiparticle traces a 
closed-loop trajectory in momentum space, its pseu- 
dospin sweeps out a certain solid angle, just as in the 
canonical Berry phase setting P2I22I Such orbits naturally 
occur when an external magnetic field B is present - 
they are the cyclotron orbits. In monolayer graphene the 
corresponding Berry phase is equal to n = \ x (2ir) at all 
energies.^ This property is the reason for the \ -shift in 
the Landau level filling factor v — Ay. [N — |) at which 



iVth magnetoresistance minimum occurs.^ES] Here the 
factor of four is the spin-valley degeneracy, assuming it 
is preserved. 

Given the unusual band structure of BLG, it is in- 
teresting to consider the effects of the Berry phase and 
other pseudospin-related phenomena on the Landau lev- 
els and t he m agnetic response in this material. Indeed, it 
is knowrfSMl that the pseudospin generates a linear cou- 
pling to the transverse component B z of the magnetic 
field, similar to a real spin. 

Note that such a pseudo-Zeeman coupling does not vi- 
olate the time reversal symmetry of the system at B = 0. 
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FIG. 1. BLG band dispersion as a function of = hvok, 
where k is the momentum measured from the nearest K ± 
point. At B = the bands are valley-degenerate. The dashed 
curves show their dispersion calculated from Eq. (JsJ) for the 
interlayer bias 2U = 240 meV. In a finite field, the bands 



acquire a pseudo-Zeeman shift, Eq. (481, opposite in the two 
valleys. The solid curves show the result for K + at B = 5T. 
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Since this symmetry operation interchanges the valleys, 
it is only the sum M+ + M~ of the corresponding mag- 
netic moments that must vanish. Further symmetry con- 
siderations require the pseudo-Zeeman shift of the energy 
eigenvalue E„ to be linear in both applied fields, 



AE q oc —E Z B Z cos30g , 



(1) 



where </>„ is the polar angle in reciprocal space relative to 
the zone center q = 0. This expression conforms to the 
following valley-interchanging operations: (i) a reflection 
O x with respect to the y-z plane, and (ii) a composite 
operation 2 consisting of a rotation through angle tt 
around x-axis in the midplane, followed by time rever- 
sal. Both of these operations leave the crystal structure 
invariant (see Fig. [3]) . The first one keeps E z the same 
but reverses the sign of B z (because B is a pseudovec- 
tor). The second changes the sign of E z but keeps B z 
the same. 

Equation ([I]) constitutes a magnetoelectric effect in bi- 
layer graphene. It implies that the valley symmetry can- 
not be broken solely by B z or by E z alone. Rather, both 
fields must be nonzero simultaneously. (It is also remi- 
niscent of the Chern-Simons term which occurs in topo- 
logical insulators.^) Below we study this kind of valley- 
symmetry breaking analytically, focusing on the question 
how it modifies the Landau level dispersion. 

Prior theoretical studies^ have already showed that 
Landau levels in bilayer graphene become valley split at 
finite U. This was explained by noting that the quasipar- 
ticle wavefunctions of the two valleys have different dipole 
moments in the z-direction. Equation (JlJ offers a com- 
plementary interpretation: the two valleys in a biased 
bilayer graphene have different magnetic moments! 34 * 35 ^ 

The ratio of the pseudo-Zeeman term ([I]) and the Zee- 
man energy due to real spin determine the effective g- 
factor of bilayer graphene. We show below that g can be 
an order of magnitude higher than its bare value g = 2. 
This resembles the situation in Bi, another low band gap 
material. In fact, there is a mathematical similarity of the 
low-energy theories^ of the two materials. (Of course, 
Bi is three-dimensional.) 

The dependence of Landau level energies in bilayer 
graphene on B z and U is known to be quite compli- 
cated (see, e.g., Refs. [38H40|. We show that it can be 
understood if one applies quasiclassical quantization to 
the Mexican hat band structure. This procedure requires 
calculating the phase shifts $ c acquired by quasiparticlcs 
on their cyclotron orbits. Both the pseudo-Zeeman term 
and the Berry phase contribute to <f> c . As a result, $ c 
generally is not an integer multiple of the monolayer value 
tt. When it does become equal to n, at certain values of 
U , an interesting phenomenon occurs: adjacent Landau 
levels of opposite valleys become degenerate. Therefore, 
there are an infinite number of Landau level crossings 
within the same band. 

Landau level crossings in the two dimensional elec- 
tron gas (2DEG) has previously attracted much theo- 
reticaP2S3 and experimental 4 ^^] interest because the 
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FIG. 2. (Color online) Graphene monolayer (left) and result- 
ing Brillouin zone (right). 



2DEG then exhibits many of the properties found in fer- 
romagnets. Therefore, BLG may be a promising system 
for studying quantum Hall ferromagnetism. 

The remainder of this article is organized as follows. A 
brief summary of BLG band structure properties is given 
in Sec. [IT} The quasiclassical approximation is discussed 
in Sec. |III| Illustrative Landau level spectra are presented 
in Sec. |IV[ The anomalous Hall conductivity of the BLG 
is computed in Sec. [V] Concluding remarks are given in 
Sec. |VI| Technical notes are gathered in the Appendix. 



II. ANALYTIC RESULTS FROM PRIOR WORK 
A. Zero magnetic field 

The band structure of BLG, well known from previous 
literature P is shown in Fig.[l] In this section we summa- 
rize its main properties, focusing on analytic results. 

The unit cell of a graphene bilayer, depicted in Fig. [3j 
consists of four atoms, which we label u, v, u, and v. 
The underlying Bravais lattice is the triangular Bravais 
lattice of either honeycomb monolayer (Fig. [2| . The Bra- 
vais lattice sites are at locations R = n 1 a 1 +n 2 <x 2 , where 

a 1 = a x and a 2 = a ^x + ^y) are primitive direct 
lattice vectors, n 1 2 are integers, and a — 2.461 A is 
again the lattice constant. The corresponding elemen- 
tary reciprocal lattice vectors are b x = 

and b 



-x 



a V3\ 2 ~ 2 J 

y. The three nearest neighbor sepa- 



ration vectors 8 l 2 3 are given by 8 1 



a 1 — \a 2l and £ 3 



3 



a, 



3 a 2' 



-g a 1 + |o 2 , each of length 
\8j\ = a /-\/3 = 1.42 A. The in-plane locations of the four 
sublattices are then given by the subscripts: u^, v-^ + g , 

^R+S ' an d ^R-S ' an d the separation between the (u, v) 
and (u,v) planes is d — 3.35 A. The (u,v) layer (B) is 
shifted by 8 1 relative to the (u, v) layer (A), a configura- 
tion known as Bernal stacking. 

Note that repeating the Bernal stacking ABABAB. . . 
generates the common form of graphite. In graphite, 
the v and u sublattices form one-dimensional chains, 
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while the u and v sites lie above and below hexagon 
centers in neighboring planes. The electronic struc- 
ture of graphite dates to the seminal work of Wal- 
lace®! and subsequent work by McClure^ and by Slon- 
czewski and Weiss p21 known as the Slonczewski-Weiss- 
McClure (SWMc) model. The SWMc model is equiv- 
alent to a seven-parameter tight binding model which 
describes nearest neighbor in-plane hopping (amplitude 
— 7o), three interplane hopping processes (71, 73, 74), two 
next-nearest plane hoppings (72,75), and an on-site en- 
ergy shift A' which distinguishes the chain sites (v, u) 
from the non-chain sites (it, v) in each unit cell. Param- 
eter A' should not be confused with A = A' + 72 — 75 . 

In BLG, 72 and 75 do not enter and one further ex- 
pecfapsl A' BLG = Ag raphitc /2. Therefore, in BLG we are 
left with five parameters: 7 = 3.0 eV, 7l — 0.41 eV, 
7 3 = 0.3 eV, 7 4 = 0.15 eV, and A' = 0.018 eV. (For 
the interpretation of these parameters within the tight- 
binding picture, see Fig. [3] For a discussion of their nu- 
merical values, including the uncertainties, see Ref. 1531 ) 
Finally, to describe a biased BLG, we include a scalar 
potential ±U on the two layers. 

The SWMc Hamiltonian of BLG in second-quantized 
notation is written as H = ^W Hq ^<j, where 



% = (4 Vq 4 4) » 



(2) 



is a four (sublattice) component creation operator with 
crystal momentum q, and 



( -U -70 S q 74 Sq 73 S *q \ 

-loSq -U + A' 7l j 4 S q 

7 4 s q 7i U + A' -7 S q 

V 7 3 S q 74 Sq ~J S* U J 



(3) 



Here, as in Ref. 160[ we define the dimensionless in-plane 
hopping amplitude 



S„ = e tqd i + e iqS * 



(4) 



In the vicinity of the two inequivalent Brillouin zone cor- 
ners q — ±K (see Fig. EjJ), S q vanishes, and writing 
q = ±K + k one finds 
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K+k 



V3 
2 

V3 



(k x -ik y )a + O(k 2 ) 



(5) 



S_ K+k = +^ Y (k x + tk y )a + O(k 2 ). (6) 

Setting all parameters but 7 to zero, one obtains the 
monolayer dispersion, 



e k = l0 \SJ=hv \k\+O(k 2 ) , 



(7) 



where v — \/3j a /2h w 1.0 x 10 8 cm/s is the Fermi 
velocity. 

If we turn on the interlayer hopping r y 1 and the inter- 
layer potential U, keeping 73 = 7 4 = A' = 0, then we 



obtain^ the spectrum 



E M 



;7i 



U 2 



A(e) 



i 7l +(7 1 2 + 4C/ 2 ) £ : 



*2 A 2 {e k ), 

1/4 



(8) 
(9) 



Here s l and s 2 label the four bands as follows: s 1 — 
± labels the conduction and valence bands, respectively, 
while s 2 = +1 for the outer bands and s 2 = — 1 for the 
inner bands. Thus, the ordering of the four levels is 



E_+ < E__ < E, < E, 



(10) 



(For aesthetic reasons, we will usually abbreviate s 12 — 
± in the subscripts, as above.) 

Due to particle-hole symmetry at 7 4 = A' = 0, we may 
restrict our attention to the conduction bands s 1 = +1. 
In this case, the shape of the energy bands is as follows. 
For the outer band, E ++ is a monotonic function of e, 

starting at q = ±K, where E ++ = E = \f^\ + U 2 and 
extending to E ++ (0) 37 (assuming 7 ^> r y 1 ,U). We 
will be interested mainly in the inner (s 2 = —1) bands, 
shaped as the Mexican hats near q = K ± , i.e., k — 
0. For example, the conduction band E + _ k has a local 
maximum — the top of the hat — at £ t = 0, where 

E^ = U and a local minimum — the bottom of the hat 

— at e k = where 



Vu 2 



E 2 , 



E^ 



AU 2 



(11) 



Hence, this minimum is attained on circles of radius k^ = 
s+/hv , centered at the zone corners. 

Inverting the relation between E and e, and suppress- 
ing the labels s 1 2 , one finds 



el = E 2 



u 2 - s 3 r 2 (E : 



k) > 



r{E) ee (75 s + AU 2 )E 2 - 7l 2 u 2 



1/4 



(12) 
(13) 
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FIG. 3. (Color online) Crystal structure of bilayer graphene. 
We label four sublattices by u, v, u, v. Also shown is the 
assignment of the hopping parameters of the tight-binding 
model. The labels ±U indicate the electrostatic potential 
energies of the layers. 
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where s 3 = ±1. This equation has no solutions when 
E 2 < E 2 , which is the band gap. There are two solutions 
when < \E\ < U, both in the inner (s 2 = —1) band. 
For U < \E\ < E = \J\J 2 + ^\ , the energy is between 
the local maximum of the inner band and the minimum 
of the outer band, and there is one solution. Finally, for 
\E\ > E there are again two solutions, one with s 2 = — 1 
and one with s 9 



-1. As we shall see in Sec. IIV Al the 



existence of two solutions E within the inner band — one 
on the inside and the other on the outside of the Mexican 
hat — gives rise to multiple level crossings when magnetic 
field is turned on. 

If one is interested only in the inner bands and low en- 
ergies, |J5 S ^1 <C Ji, then one can implement a unitary 
transformation, discussed in Appendix |Aj which decou- 
ples the inner (s 2 = —1) and outer (s 2 = +1) bandsPEH 
The results of this procedure are further described in 
Sec. HTB21 
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B. Quantizing magnetic field 

In the presence of a magnetic field B = Bz, the two 
components of the wavevector no longer commute. Near 
the zone corners, we invoke the Kohn-Luttinger sub- 
stitution fcj — > TTi such that [7r,,.,7r„] = —i/£ B , where 

l B = \Jhcje\B\ is the magnetic length. We define the 
ladder operators, 



(n x — in y ), 



ot = - 



(n x + iiTy) , (14) 



which satisfy the commutation relation 

[a, at] = sgn(B). (15) 

Using Eqs. ^ and @, wc find the Hamiltonian of K + 
valley to be 



H + (U) 



77 4 w a 



— w a 
■U + A' 
7i 

t 



7i 
U + A' 



774 uj q a 
—uj n a 



u i 
U 



J 



where ?y 3 = 7 3 /7 = 0.1, tj 4 = 7 4 /7 = 0.05, and 



35meV v / |-B(T)| 



(16) 



(17) 



Throughout we shall ignore the effects of real Zeeman 
splitting, which are small due to the value of the Bohr 
magneton, n B = eh/2m c c = 57.9/xeV/T. At the highest 
fields in the relevant experiments (B m 30 T) the real 
Zeeman splitting is on the order of a few millivolts, which 
is much smaller than even the smallest of the SWMc 
energy scales. (As we shall see, the pseudo-Zeeman effect 
can be significantly larger.) 

The Hamiltonian H~ of K~ valley is obtained from 
H + via the replacements 



Rx 



(18) 



FIG. 4. (Color online) Landau level energies vs. magnetic 
field for U = (left) and U = 80meV (right). Solid lines cor- 
respond to the K + valley and broken lines to the K~ valley. 
The color distinguishes the spectra of "H a (black), "Hb (red), 
and H c (blue), where H a ,b,c are defined in Appendix B. 



which is the reflection in the y-z plane. The commu- 
tation relation (15) between a and a* and therefore the 



energy spectrum is preserved if we additionally reverse 
the magnetic field, 



Rb 



B 



-B. 



(19) 



Taken together, these replacements implement the sym- 
metry operation 0% = R x Rb discussed in Sec. [TJ The 
other valley-interchanging operator 2 is represented in 
terms of the unitary matrix 



V 



a x 

7 X 



(20) 



and the time-reversal operation Sq — > (Sq)*, i.e., 

R T : a-t-a, a f -> -a) , B -B . (21) 



It is easy to see that 

H~{U) = R T R B [V f H + (-U) V] 



(22) 



Since RtRb also does not change the commutation rela- 
tion (15), the spectra of T-L~{U) and H + (—U) coincide. 
Thus, it suffices to discuss the spectrum of T~L + , from 
which one can obtain the spectrum of TAT by reversing 
the sign of either U or B. 

These symmetries further imply that at B = the two 
valleys are degenerate in energy and that additionally, 
each valley is symmetric under U —>—[/. On the other 
hand, at finite B, the valleys are degenerate only if U = 0. 
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Note also that the total spectrum, including both valleys, 
is particle-hole symmetric when A' = 7 4 = 0. 

Making use of the eigenvectors \n) of the number op- 
erator a^a, we write the general bilayer wavefunction as 




!*> = £ 




(23) 



The matrix representation of the corresponding Hamilto- 
nian is discussed in Appendix B. If all SWMc parameters 
are kept, it can be diagonalized only numerically. Some 
results are shown in Figs. |3J [5] and [6j which illustrate 
that the spectrum can be rather complicated. In the re- 
mainder of this section we review certain limits where 
some analytical progress can also be made, which helps 
with physical understanding of these results. 




40 80 
U (meV) 



120 



1. 73 = limit 

It is simplest to consider the case where 7 3 = 0, which 
turns out to be an excellent approximation at large fields. 
When 7 3 = 0, the eigenstates of W + fall into one of three 
classes: 



FIG. 6. (Color online) Landau level energies vs. interlayer 
bias U for a field value B — 5 T. Solid lines correspond to 
the K + valley and broken lines to the K~ valley. The color 
distinguishes the spectra of "H a (black), Hb (red), and H c 
(blue). At an accidental degeneracy (level crossing), the color 
and the line type cannot both be identical. The shaded area 
indicates the energy gap at B = 0. 



\A) = 




(24) 




FIG. 5. (Color online) Landau level energies vs. magnetic 
field B for the case 73 = 74 = U — 0. The labeling of the 
states corresponds to that in the text. 



and 



/u n _ x \n - 1) 

u n \n) 
\v n +i \ n + !) 



(25) 



with n > 1. Clearly jV'-i) is an eigenstate with eigen- 
value E = U. Applying H + to \ip ), one obtains the 3x3 
Hamiltonian for the ip sector, 




r) 4 uj 



U 



(26) 



Finally, the spectrum for the \ip n ) states (n > 1) is given 
by the eigenvalues of the 4x4 Hamiltonian 



H„ = 



where 



—U 


—W 







-w n 


-17 + A' 


7i 






7i 


U + A' 


-w n+1 







-w n+1 


u 




W = 







(27) 



(28) 



To label the states, it is helpful to consider the case 
7 3 = 7 4 = A' = U = 0, corresponding to a pure nearest- 
neighbor hopping model with constant (zero) local site 



6 



energies. One then finds the following (valley-degenerate) 
spectrum: 




-n, Sl ,s 2 =<VM/o + l^+o ) P + s 2Cn > (30) 



where 



and 



7i 



136 T 




(31) 



(32) 



Note that the n = — 1 and n = states require a separate 
labeling convention. 

For the full model, with 7 3 , 74, and A' restored, 
particle-hole symmetry is broken by the 74 and A' terms. 
These are relatively small however, so there remains an 
approximate particle-hole symmetry, as shown in Fig. [6] 
The state labels are then defined by adiabatic continuity 
with the 7 3 = 7 4 — A' = limit. 



2. Low energy effective theory 

As mentioned above, at low energies one can imple- 
ment a unitary transformation, which decouples the in- 
ner (s 2 = — 1) and outer (s 2 = +1) bands or der b y or- 
der in S, which vanishes at the zone corners Here 
S = =f(v / 3/2) ao(k x =p iky) and the upper (lower) sign 
denotes the K + (K~) valley. To order S 2 one obtains 



H = 



\{A + 2U) SS+ -V 



7 3 



,S't 



7i 



\ 



V 



73 



s _ To 5 t2 



7i 



X(A-2U) S^S + U 



where A = (7o/Ti) 5 



A = A' 



53.5 and 

, 2 7i7 4 
7o 



59 meV 



(33) 



(34) 



is a composite parameter describing electron-hole sym- 
metry breaking effects of A' and 7 4 . Anticipating the 
introduction of an external magnetic field, we have al- 
lowed for the possibility that S and S> do not commute, 
cf. Appendix [A| for derivation. 

The eigenvalues of H to leading order in 7 3 are^ 



,-,h 



f 2 ~ 
%A 
Ti 




^icos3^, 
7o7i 



(35) 



where <p is the polar angle of fc. This agrees with Eq. 
in the appropriate limit. 



The 2x2 form of matrix H in Eq. ( 33 1 naturally leads 



to the concept of a pseudospin- ^ degree of freedom which 
simplifies calculations somewhat. We use this approach 
sparingly for the following reasons. First, in experiments 
U is not necessarily much smaller than r ) 1 , in which case 
the reduction to a two-band effective Hamiltonian is not 
valid. Second, the calculation of the pseudospin-related 
effects are not difficult even when all four bands are kept. 
Finally, the low-energy theory does not produce accurate 
results for the Berry phase. A brief discussion of this 
technical issue is also given in Appendix [A} 
In a nonzero magnetic field, H + becomes 



( B(A + 2U)aa) ~-U ^ w at - /3 7l a 2 \ 
7o 



1 ^w a-/3 7l at 2 0(A - 2U) da + U 

\ Tin 1 



(36) 



To 



while H~ is obtained via substitutions (18 1. When 73 = 
0, their eigenvalues are easily obtained by considering the 
basis of states 



+ \ - 



+ 

n 

+ 'n+i; 



n 



(37) 



In this basis the above Hamiltonian takes the form 
f fi (A + 2U) n-U -P^n{n + l) 

-/3 7lV /n(n + l) p (A - 2U)(n + l) + U, 



(38) 



When ri = —1, we have = and the energy levels 
in the two valleys are £_j = ±Z7. With n = 0we again 
havc Uq =0, and 



= f3A ± (1 + 2/3)17 . 



(39) 



The splitting of the n = —1 and n = levels and 
their valley-dependent slope as a function of U lead to 
a characteristic diamond-shaped crossing pattern, shown 
in Fig. [7] The largest energy gap occurs in the unbiased 
sample, U = 0, and its magnitude s» 0.5 meV x B(T) is 
comparable to that measured in Ref. 27 in a suspended 
BLG. On the other hand, an order of magnitude smaller 
gaps (smaller than even the bare Zeeman gap) have been 
observed in a more disordered BLG on Si02 substrate!^! 
Finally, for n > one has (similar to Ref. |40|) 



± 



1 



r BA i 2 

(2n + l)3U T ~ U\ + n(n + 1) B 2 ^ 2 



(40) 

This completes our summary of the (mostly) known an- 
alytic results for the energy spectrum of BLG. 
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III. QUASICLASSICAL APPROXIMATION 
A. Effective g-factor 

Renormalization of the electron magnetic moment is a 
well-known phenomenon in the solid-state physics. Most 
often it comes from spin-orbit interaction; however, in 
crystals without inversion symmetry there is an addi- 
tional contribution due to the orbital angular momen- 
tum: 



(a\M\a) 



e 

2c 



air x v \a) 



(41) 



Here a is a given Bloch state and r, v are the position 
and velocity operators, respectively. Since we are not in- 
terested in the center-of-mass motion, in evaluating M a 
we must assume that the expectation value of position 
vanishes, i.e., that r has only off-diagonal matrix ele- 
ments^ 3 



a ^ a' . 



This leads to 



M « = ^E[HV fe |c/)x(«>|a)] 



(42) 



(43) 



A lucid derivation of Eq. ( 43 1 was given previously in 
Refs. [55] and which also contain references to much 
earlier workP^l 

Below we assume that B and M are both in the z- 
direction. The orbital contribution to the <?-factor is 



zs: 









N N 
\ \ 




Sn = 


\ \ 

N / ' 




y/n = -1 












\n = -1 




V 








s X n = o . 

\ 

i , \ 



-4 -2 2 4 

U (meV) 

FIG. 7. (Color online) Landau levels —1 and as a function 
of U at B = 5T, i.e., f3 = 0.037. The solid (dashed) lines 
correspond to K + (K~) valley. The SWMc parameters are 
taken from Ref. Eland A = 59meV, cf. Eq. p4l. 



g = 2M a /fi B where fi B — eh/(2m c c) is the Bohr magne- 
ton and m e is the bare electron mass. To calculate M a , 
we can add and subtract the omitted diagonal term in 
Eq. ( 43 ) , which gives 



where 



F a = -i(a\V k x v\a) ■ 
D a = -i [(a\ V fc | a) x v 



(44) 



(45) 
(46) 



(note that both F a and D a are real) and where 

H-^Ea (47) 



(a v \a) 



is the group velocity vector (the subscript a in v g is omit- 
ted for simplicity) . Using these formulas we compute the 
energy dispersion 



E n - BM n 



(48) 



It is interesting to compare our formula with those in 
literature. A very close analogy is provided by Bi, whose 
effective Hamiltonian is also a 4 x 4 matrix linear in k. In 
an early paper^ 3 where the calculation of the g-factor of 
Bi is discussed, the subtraction of the diagonal term D a 
is lacking, so that the result is not gauge-invariant. Below 
we show that D a is related to the Berry phase, which ap- 
parently has not been handled correctly in Ref. [37| (con- 
sidering that it precedes Berry's worlpSJ by almost two 
decades, it is hardly surprising). 

Let us now apply our general formula to BLG. For K + 
valley we can choose the eigenvectors of H + in the form 



) = (u a e~ 



I'o 



U 



itp\T 



(49) 



where a now throughout this section stands for 
{s 1 ,s 2 ,fc}. It is assumed that the imaginary parts and 
the entire dependence on <p — the polar angle of k - 
enter via the exponential factors only. A straightforward 
calculation yields: 



F, 



+ _ 



2v n 



2kE a 



4e 2 



D 



+ _ 



2s 2 A^e k ) / 
2hvlU 
M 2 (^) 



(50) 



(51) 



where A(e) is given by Eq. flob . The eigenvectors for K + 
valley can be obtained by replacing e ±iip in Eq. (49) with 



-e Ttv and so the signs of F a and D a are reversed. 
The last term represents the pseudo-Zeeman effect due 
to the orbital magnetic moment. Algebraic manipula- 
tions with Eqs. p4| , ( |5"0| ), and ([51]), together with 
the relations 



1 dE 
h dk 



£ S- 



■,r 2 (E) 

E s 2 A*(e) 



(52) 




0.2 0.3 0.4 
£fc (eV) 



superscript denotes the valley, as usual): 



at 



Ml 



= ± 



8m e v 2 U 



t'B 



7?" 
14- 



4t/ 2 



(56) 



22. This 



Thus, at U — 100 meV we obtain 
is one order of magnitude higher than the bare value 
g = 2 and is about as large^ as in Bi. (For this rea- 
son, we neglect the bare Zeeman coupling in this arti- 
cle.) For U -C 7i, the effective g-factor is proportional 
to U, as appropriate for the linear magnetoelectric cou- 
pling [Eq. Q]. Therefore, a roughly linear variation of 
the band edge positions with B and U can be expected. 
This issue is addressed in more detail in Sec. IIVI 



FIG. 8. (Color online) Orbital magnetization M+_ of K + 
valley as a function of Ek = hvok at U = 0.1 eV. The location 
of the band bottom e k = e+ is marked by the arrow. 



and 



yields 



Ml 



s 3 r 2 (E a ) - s 2 A 2 {e k ) = -7? +2/7 2 , 



eh 



c 7 f + 4( 7 2+4f/2) £ 2 



El 



(53) 



(54) 



For the lower energy conduction band, on which we 
mostly focus later, M+_ k is plotted in Fig. [SJ The mod- 
ified spectrum E a is plotted alongside E a in Fig. [l]for all 
four bands and in Fig. [9] for the lower conduction band 
only. At k = we have a particularly simple result, 



9 Sl ,s 2 ,o 



M, 



Mb 



2 U 

o 7 2 



(55) 



for all s 1 2 , in agreement with Eq. (54) of Ref. [35] 

As one can see from Fig. [8j the ^-factor has an in- 
triguing energy dependence, which prompts the question 
of whether it can be verified experimentally. Unfortu- 
nately, this appears problematic. There is no optical 
transition between the energy levels split by the pseudo- 
Zeeman effect as they belong to different valleys, and so, 
methods analogous to the electron spin resonance would 
not work. Another conventional method of extracting 
the (/-factor would be to measure the valley-splitting of 
the Shubnikov-de Haas effect. However, this splitting 
also includes the contribution of the Berry phase, dis- 
cussed later in this Section. This contribution effectively 
compensates for nonmonotonic variation of the g-factor, 
making the valley-splitting of Landau levels only weakly 
dependent on the Fermi energy (or Landau level index). 

The most easily observable manifestation of the 
pscudo-Zceman effect appears to be the displacement of 
the band edges, e.g., the bottom of the Mexican hat of 



the conduction band. At this point, Eq. (54) yields (the 



B. Quantization rules 

While numerical calculations of the Landau level spec- 
trum is possible for any choice of parameters, in Sec. [IV] 
we shall see that the result can be rather complicated. 
Therefore, both exact and approximate analytical meth- 
ods remain valuable for this task in hand. So far, we have 
discussed two such methods. First, for U = 7 4 = A' = 0, 
closed-form expressions for the Landau level energies 
[Eq. ([9])] exist. Second, if these energies are much smaller 



than 7 1; then the approximate Eq. (401, valid for finite U, 
can be used. In this section we outline another approach 
— the quasiclassical quantization — which can be used 
for arbitrary relation between U and 7l . Within this ap- 
proximation, Landau level energies £^ s are taken to 



be equal to the renormalized band energies ( 48 ) evalu- 



ated at certain quantized orbits in the reciprocal space: 



n.s-. ,s q 



(57) 



If we ignore 73 , the orbits are circular and the area of the 
nth such orbit satisfies the Onsager condition 71 



7r(**4) a =27r(n + ^) 



(58) 



where k^ is the radius of the orbit and d„ is a dimen- 
sionless number discussed below. 

The quasiclassical approximation is accurate through 
the order 0(£g 2 ) or alternatively 0{l/n). It turns out 
to be exact for parabolic dispersion (where S = i) and 
in monolayer graphene (where 8 = 0). The quasiclassi- 
cal approximation for general matrix Hamiltonians was 
previously studied in Refs. [721 and 1731 and specifically in 
the context of graphene in Ref. [TJl However, we found it 
most instructive to follow Refs. l68l and l75l 

The physical picture is as follows. In a weak magnetic 
field, momentum k of a quasiparticle slowly rotates as a 
function of time t according to the equation of motion 

2-7T 

k = u c z x k , w c = — sgn(w g ) , (59) 

where T = 2irk n i] 3 /\v g (k n )\ is the cyclotron period. (For 
simplicity, the valley and band labels are temporarily 
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omitted.) The rotation of fc causes a slow evolution of the 
wavefunction \a) in the pseudospin, i.e., sublatti ce spa ce. 
This causes the accumulation of the Berry phase^ES 



$ B = sgn(> ) / dt (a \iV k \a) ■ k . 



(60) 



Thus, in monolayer graphene where $ B = n, we get S = 
0, which implies the existence of a level at zero energy.^ 
Comparing Eqs. (46) and (60) we see that for the 
isotropic spectrum, i.e., for 7 3 = 0, we have 



2nk 



(63) 



The quasiclassical quantization rule i: 



68 



sgn(w K ) <p dny t B tt x + ^ b = n(kJ B Y + $ B 

= (2n+ 1)tt. 



(61) 



This formula can be understood as a generalized Bohr- 
Sommerfeld rule: since £ B ir x plays the role of "momen- 
tum" conjugate to the "coordinate" Tr y , the top line rep- 
resents the total phase shift acquired along the orbit, 
including the geometric phase. Equation (61) establishes 
the precise relation between the Onsager number 5 and 
the Berry phase <I> B : 



S = 



1 



2tt ' 



(62) 




FIG. 9. (Color online) Evolution of a particular (n = 5) 
Landau level of the K + valley as a function of U. Super- 
imposed are the spectra at zero field (thin traces) and that 
with pseudo-Zeeman correction in a magnetic field B — 5 T 
(thick trace), (a) At small U, the quantized cyclotron orbit 
is outside the Mexican hat. (b) For certain U, the orbit goes 
inside the gap of the zero field spectrum, (c) At larger U, it 
moves underneath the Mexican hat where the direction of the 
group velocity is opposite to the momentum, (d) At very large 
U (not presently accessible in experiments), where the BLG 
spectrum consists of two copies of monolayer spectra shifted 
by ±17, the 71 th electron Landau level of BLG approaches the 
(n + l) th hole Landau level of the higher energy monolayer. 



Postponing the discussion of this equation for just a mo- 
ment we note that for v g 7^ 0, another version of the 
quantization rule can be established.^ To this end one 
defines a modified orbit radius k n such that 



f 



E 



(64) 



To the leading order in £>, the rule that determines k n is 
similar to Eq. (61) except 3> B is replaced by a different 
phase shift <& c : 



ir(k n l B ) 2 = (2n+l)^-$ c 



MB 



irk 



T=—(F a + D c 



(65) 
(66) 



With further analysis it is possible to show that our 
$ c coincides with the "semiclassical phase" defined in 
Ref. EH Therefore, the difference between $ c and $ B 
noted in that paper is entirely due to the pseudo-Zeeman 
shift rather than a violation of adiabaticity. 

Applying the above formulas to BLG, we obtain 



B 

2tt 



U 
2~K 



44-7? 



V7f+ 4(7? +4^)4 



(67) 



At finite U the Berry phase is a nonmonotonic function of 
momentum, which is addressed in more detail in Sec. [V] 
Here we comment only on the simple case U — 0, where 
Eq. (67) gives <1> B = at all k ^= 0. This seems to contra- 

±2ir made in most of the 



diet to the assignment 



previous work! 2 ! 13 ! In fact, there is no contradiction be- 
cause the Berry phase is not unique: different choices for 
an overall phase of the wavefunction in Eq. ( |60| can shift 
$ B by an arbitrary integer multiple of 2-7T. In the context 
of Landau quantization, such shifts can be compensated 
by relabeling the Landau index n, so that the physical 
quanitities — the radii k^ of the orbits and their energies 
— remain the same. 



Combining Eqs. (66) and (67), we obtain the analytic 



formula for the semiclassical phase: 



2tt 



Ts 3 



UE a 



V(7i 2 +4C/ 2 )i? 2 ~7? C/ 2 



(68) 



This equation should be used away from momentum k+ 
where its denominator vanishes. Finally, the quantiza- 



tion rule ( 65 ) becomes 



2tt 



(69) 
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the precise relation between n and £ a for 
'., = A' = reads 



In comparison, 

the case 7 3 = 7 4 = A' = reads 



£ I + u 2 



2ue„ 



(70) 



This result follows from Eq. ( 30 ) ; the composite label a 



denotes the set {n, s x ,s 2 }. The semiclassical Eq. (69) 



does agree with the exact Eq. (70) to the leading order 



in Wg, i.e., 0(l/n) at large n. Fortuitously, it is also 

valid for n = — 1. It predicts k^ = 0, which entails 

£i__ = U, in agreement with our earlier result. 

The valley splitting of the Landau levels can be ex- 
pressed as follows: 



-£- 



2h 



SlS 2 



^ 1 l + A{ 1 (+AU 2 )el 



(71) 



Here either fc„ or fc„ can be used in place of k because 
this formula is valid only to the leading order in (3. At 
low energies, it simplifies to 



c+ _ c- 



-2/3Z7, n>l. 



(72) 



in agreement with Eq. (40). We see that unlike the 



pseudo-Zeeman term, discussed in Sec. |III A[ the net 
valley-splitting of the Landau levels has little energy or 
n dependence. 

It is now straightforward to apply the above quantiza- 
tion rules in order to understand qualitatively the evolu- 
tion of some n ^ 1 Landau level as a function of U. For 
the K + valley, illustrated in Fig. [9j the situation is as 
follows. As U increases starting from zero, the radius of 
the orbit changes only slightly because <£> B /27r ~ 1 <C n 
for all U . On the other hand, the Mexican hat expands 
in both height (energy) and width (momentum). As a re- 
sult, the quantized orbit slips from the exterior (k > k^) 
to the interior (k < k+) of the hat. In the process, the 
orbit passes through a region where its energy is inside 
the gap of the B = dispersion because of the negative 
pseudo-Zeeman term. [For the K + valley this occurs only 
if U > but not if U < 0, see Eq. ([56) and Sec. [TV] be- 
low.] Eventually, at very large U, the orbit approaches 
the n + 1st hole Landau level of graphene monolayer, 
except it is shifted upward by U. 



IV. LANDAU LEVEL SPECTRUM 

A. Level crossings 

In this section we explain the physical origin of a non- 
monotonic [/-dependence of Landau level energies, which 
gives rise to a complicated net-like pattern with numer- 
ous crossings, see Figs. [T0| and |T5) It should be clarified 
that electron interactions, which are ignored in our cal- 
culations, can produce significant corrections to the Lan- 
dau level spectrum. However, we expect that topological 
properties of the level diagram would not change much. 



Figure 10 shows the first several Landau levels, which 
we calculated numerically as a function of U at a rep- 
resentative magnetic field of B = 5T. Only U > are 
shown because the energies at negative U can be obtained 
from the symmetry relation £+ (U) = £~(—U). Let us 
focus on the s x = +1 levels and consider the limits of 
small and large U (a similar argument can be applied to 
the s 1 = — 1 levels with appropriate sign changes). 

For small U, the Landau levels are roughly equidistant 
and those with higher index n have higher energies (in 
agreement with Eq. (30) for U = 0). In the opposite 



limit of U 3> 7i, from Eq. (16) it is easy to see that the 
BLG spectrum consists of two copies of the monolayer 
spectrum shifted by ±U. Accordingly, the set {£ n+ _} 
approaches the Landau level energies of the holes in the 
monolayer P but shifted by U > 0: 

£n+- — U — y/n~+l ui , £,7+- — U - Vn oj a . (73) 

In this limit states of higher index have lower energies. 
Therefore, any two levels of the s 2 = — 1 band cross at 
some value of U. This occurs when the corresponding 
quantized orbits are located at the same energy but on 
the opposite sides of the Mexican hat (see Fig. [9| . 

In addition, it is possible to have crossings of orbits 
on the same side of the Mexican hat if they belong to 
opposite valleys. In the semiclassical approximation, this 
occurs whenever "^/V is an integer. In this case the 
difference in n is compensated by the difference in the 
semiclassical phase, yielding the same momentum k n and 
energy E Si ^ 2 (k n ) (see Eq. (65)). For example, at U = 
we have $^ = 0, so that all Landau levels should be 
(and are) valley-degenerate. Next, |$^| — > n as U — > oo, 
so in this limit the adjacent Landau levels coincide, in 



agreement with Eq. (73). Using Eqs. (11), ( 13 ), and (68 1, 
one can show that the condition |<f>! t | = Nir is met at 



£ 2 - 



El 



1 - (2EjN 7i y 



( 73 = 7 4 = A' = 0) . (74) 



This implies that the level crossings are confined to the 
range of energies E+ < \£ a \ < \U\, which is precisely the 
range between the top and the bottom of the Mexican 
hat. The crossings at the top of the hat are between 
the adjacent Landau levels (N = 1). Since the special 
level £ d: 1 = ±U also happens to be at the same energy, 
these are actually triple crossings. In the Sj = ±1 band, 
they involve nth level of K ± , the n — 1st level of K T , 
and the —1 level of K (assuming U > 0). When j 3 = 
7 4 = A' = these unusual triple crossings appear when 
U = U n = h \/noj Q . We can show by algebraic means that 
finite 7 4 and A' give corrections to U n but do not lift the 
triple degeneracies. We suspect that this property stems 
from some hidden symmetry of the Hamiltonians . 



B. Trigonal warping 

The parameter 7 3 has a number of interesting effects 
on both the zero-field and Landau level spectra. It mixes 
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Landau levels of the same valley with indices n different 
by an integer multiple of three, see Appendix [B] This 
turns crossings between such levels into avoided cross- 
ings. Strictly speaking, we can no longer label Landau 
levels by {n, s 1 , s 2 }. Nevertheless, the effect of 7 3 is small 
at low U, so that with proper care it is possible to track 
the levels through the avoided crossings and still retain 
our labeling scheme. The calculation of the Landau level 
spectra with 7 3 ^ is handled numerically. To account 
for the level mixing at high U we had to diagonalize ma- 
trices of size 4J with large enough J ( J s» 100) to ensure 
numerical accuracy, see Appendix [B] One effect of 7 3 is 
to lift the triple degeneracy of the adjacent Landau levels 
by moving the crossing point energy away from the top 
of the Mexican hat, as expected. 

A more interesting effect is the shift of the B = band 
edges, which are the boundaries of the central band gap 
m Fig.U This can be understood as follows. The hop- 
ping 7 3 induces a trigonal warping of the zero-field bands, 
as described by Eq. (35). Accordingly, the low-energy 



region of the conduction band develops three kidney- 
shaped pockets along k = k^ circle centered, in K + val- 
ley, at (p = ~tt, 7r, and |tt angular positions.^ To the 
leading order in 7 3 , their energy is lowered below E+ by 



5E 



V8J3 U 2 
7o 7i 



7 3 U 

— < — < 1 

7o 7i 



(75) 



which follows from Eqs. (|35j). Accordingly, the band edge 
of the conduction/ valence band at B = shifts by ^fSE. 
For example, at U — 0.15 eV we obtain SE « 8meV. 
This is in a good agreement the numerical results shown 
in Figs. [6] and 10 



The effect of 7 3 on Landau levels is even more striking. 
As one can see from Fig.[l0j it leads to a bunching of Lan- 
dau levels near the conduction (and valence) band edges 
as U increases above 0.1 eV. Apparently, these Landau 
levels, which can be labeled n = — 1, n+, and % + 1, 
become nearly degenerate. Within a simple quasiclassi- 
cal picture, the explanation is straightforward: this trio 
of levels correspond to three orbits, which are identical in 
shape and energy but are separately confined inside the 
three equivalent pockets. 2 In a more refined description, 
such orbits are hybridized by a weak quantum tunnel- 
ing, so that the Bloch functions have equal amplitude in 
each pocket but different phases. To verify this picture, 
we chose a set of U in the range between and 0.15 eV 
and for each of them computed the Bloch function of 
the lowest-energy state numerically. We took 7 3 = 0.15, 
for which there is only a single threefold degenerate level 
lying just within the central gap. At all U, these func- 
tions exhibit maxima centered at ip = ^ir, it, and |7r, 



as expected (see Fig. 11). However, for U>0.1eV such 
maxima become very sharp, consistent with the picture 
of confinement and in concert with the coalescence of the 
energy levels into a single narrow bunch, as in Fig. 
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In general, the influence of 7 3 on the spectrum gets 
stronger as B decreases or U increases. This is because 
the depth SE of the pockets and their width increases 
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FIG. 10. (Color online) Landau level energies vs. interlayer 
bias U for a field value B — 5 T. (a) Top panel: 73 = 0; (b) 
bottom panel: 73 = 0.3 eV. The color and line type are as in 
Fig. [6] Note the bunching of levels at the edges of the central 
band gap when 73 7^ 0: the two levels just below the gap for 
t/;> 100 meV are both very nearly threefold degenerate. 



with U while the area in momentum spa ce p er orbit is 
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Hence, at 



equal to 2ir/£g oc B, as discussed in Sec. 
large U and/or small B, each pocket may host several 
orbits, so that higher-energy Landau levels can also form 
bunches of three, as is apparent in Fig. [lOj where there 
are two nearly three-fold degenerate sets of Landau levels 
separated by 10-20 meV from a tangle of higher energy 
states. The first bunch emerges at U ~ 80 meV and the 
second at U « 120 meV. Conversely, as B increases at 
fixed U, separate orbits no longer fit into the pockets 
and they unite into a single contiguous loop. At this 
point, the effect of 7 3 can safely be neglected. 



C. Energy gap 

The above discussion indicates that the energy gap of 
BLG can be controlled not only by U but also by B while 
keeping U fixed. Since this gap can strongly affect the 
low-temperature transport, it may be of interest in appli- 
cations, and so it deserves some discussion. The magnetic 
field tends to reduce the gap relative to the zero-field 
case, as one can see in Figs. [6] and 
area indicates the zero field gap. In other words, 



where the gray 
some 

Landau levels can reside inside the bandgap \E\ < E^ of 
the B = spectrum. This phenomenon is a direct man- 
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ifestation of the pseudo-Zeeman shift. It is seen more 
clearly in Fig. |9][b), where only one Landau level (from 
the K + valley) is shown. For a certain U > this level 
drops below the zero-field minimum of the conduc- 
tion band. Similarly there is another Landau level from 
K + valley not shown in the Fig. j9^b), which rises above 
the maximum — E+ of the valence band. This is because 
the pseudo-Zeeman effect has opposite signs in the two 



0.16 



valleys. Based on this argument, we can use Eq. (56) to 
show that, e.g., the bottom of the conduction banc 
to 



shifts 



1 + {2Uh 1 ? 



(76) 



where (3 is defined in Eq. (31 ). In principle, this approx- 



imate formula can be refined by semiclassical quantiza- 
tion. The true band edge is determined by the lowest- 
energy Landau level of the conduction band. Its index 
n^, which depends on U and B, can be found by setting 
k = k^ and dropping the second term on the right-hand 
side of Eq. m-. + § ~ ej/w§ = (l//?)(ej 7 i) 2 - A 



similar result can be obtained from the low-energy effec- 
tive theory, by minimizing the energy in Eq. ( |40[ ) with 
respect to the Landau index n. With A' = we obtain 



1 



2U 2 



13 7l 2 + 4C/2 



(77) 



Since (3 cx B, our approximate formula 2Ef for the gap 
predicts a linear gap narrowing as B increases at U = 
const. Figure [12] demonstrates that it is quite accurate 
up to the point where nj drops to zero, i.e., up to the 




FIG. 11. (Color online) Absolute value of the Bloch function 
for the lowest-energy Landau level of the conduction band. 
The origin is at the K + point, the radial coordinate is kin, 
and U= 0.15 eV. 




FIG. 12. (Color online) Energy gap separating Landau levels 
of the valence bands from those of the conduction band as a 
function of the magnetic field B. The cusps on the curves are 
due to discrete changes in Landau level index n* (see main 
text). The upper solid curve is for 73 = and the lower one 
for 73 = 0.30 eV. The analytic estimate per Eqs. (76 1 and 
(781 is shown by the dashed line. 



field where j3 « 2U 2 /j 2 . Of course, this approximation 
misses the small cusps produced by the discrete changes 
in ?v 

At larger B, the gap is determined by the energy of 
the special n = Landau level for which Eq. ( 76 ) is not 
valid. Instead, we can use Eq. (39) to get 



E+ - (2\U\ - A) (3 , /3<2C/ 2 /7i 2 



(78) 



We see that the -B-dependence remains linear but the 
slope becomes larger by a factor of two or so. This pre- 
diction is in a reasonable agreement with numerical cal- 
culations (Fig. 12). The deviations seen at B > 10 T are 



due to insufficient accuracy of the low-energy theory at 
such fields. The total reduction of the gap as the field 
changes from B = 0T to 15 T is about 15meV or 10%. 

At even larger B, level n — on the s 1 = sgn(f) side 
would cross with level n = — 1, so that the slope of the 
linear dependence would change again. That n = — 1 
level would eventually intersect with the other n = 
level if B keeps increasing, at which point the gap would 
momentarily vanish. An example of such an intersection 
is shown in Fig. [7] (although the energies are plotted as 
a function of U). 

Let us now discuss the effect of 73. In Fig. [6j the ener- 
gies of the lowest-energy levels of the conduction and its 
counterpart in the valence band seem to be lined up with 
the respective edges of the B = spectrum, as though 
the pseudo-Zeeman effect is canceled. This cancellation 
is fortuitous. We attribute it to the zero point motion 
of the orbits confined inside the pockets. Clearly the 
Bloch functions (Fig. |11[ ) have some finite spread around 
the centers of the pockets. Thus, in the conduction band 
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such orbits are raised in energy above the actual minima 
of the band, which counteracts the effect of the pseudo- 
Zeeman shift. Indeed, a better measure of the pseudo- 
Zeeman effect is the valley splitting, which is nearly the 
same in Figs.flO^a) and (b). (The latter is essentially the 
upper half of Fig. [6). The magnitude of the zero-point 
energy shift depends on U and B and just happens to be 
numerically close to the pscudo-Zeeman shift in a range 
of parameters used in Fig. [6j 

The gap narrowing becomes more apparent at larger 
magnetic fields, see Fig. 12. The upper and the lower 
solid curves represent the energy gap without and with 
7 3 , respectively. At B = 0, the distance between the two 
curves is about 8 mcV, which is close to 2 SE « 9 meV 
per Eq. (75 1. As B increases, this distance quickly di- 



minishes, and the gap concomitantly narrows. 



V. ANOMALOUS HALL CONDUCTIVITY AND 
VALLEY MAGNETIZATION 



Systems that carry a finite Berry phase can exhibit a 
nonzero Hall conductivity <7 H even at B = 0. This is 
known as the anomalous Hall effect (AHE). The AHE 
and other manifestations of the Berry phase in electronic 
properties have been recently reviewed in Ref . 1771 It has 
been shown that for a partially filled band, <r H in units 
of e 2 /h is equal to the Berry curvature 



x la i 



IV, 



(79) 



integrated over all occupied states. By the Stokes' theo- 
rem, in the two-dimensional case the result is determined 
solely by the Berry phase at the Fermi level. Therefore, 
we can readily compute the anomalous contribution to <j h 
from our Eq. (67). To do so, we need the Berry phase as 



a function of energy. Substituting Eq. (53) into Eq. (67), 
we obtain 



<I> 



B 

2tt 



U 2E 2 



7l 2 -2 S3 r 2 (£) 



E AU 2 



•7? 



2s 3 r 2 (E) 



(80) 



where r(E) is defined in Eq. (13). The opposite signs in 
this formula indicate that the two valleys give opposite 
contributions to the AHE. Therefore, tr H is nonzero only 
if unequal population of the valleys is created. While 
this occurs naturally for B ^ 0, we desire, in the con- 
text of the AHE, that it should also occur in the absence 
of an external magnetic field. Theoretical proposals for 
achieving that have been advanced in Refs. |34l I78H801 
Here we do not address any mechanisms of valley po- 
larization but simply compute all the quantities for K + 
valley only. Comparison with previous work will be given 
at the end of this section. 

For brevity, we limit the consideration to the case 
where the Fermi level [i resides in the conduction bands 



{s\ = +), i.e., /i > 0. Using Eq. (67), we obtain: 



where 





<Mm) - Mm) 



E,< f i<U, 



a-{n) - &+(U) + ct+(m) E <fi, 



a S3 (E)=^$+(s s ), 



(81) 



g s = 2 is the spin degeneracy, and 



£ o ^£ ++ (fc = 0W 7l 2 + £/ 2 



(82) 



(83) 



Hence, E^, U, and E are the energies of the Mexican 
hat bottom, Mexican hat top, and the upper conduction 
band bottom, respectively. At these energies the topol- 
ogy of the Fermi surface changes: from two concentric 
circles to one and back to two (we ignore 7 3 ). Accord- 
ingly, CT H (/i, U) is nonanalytic at such energies: it has 
discontinuous derivative (cusps), which are marked by 
the dots in Fig. [l3^a). Note that in the limit of small 
U, er H approaches the universal value of 2e 2 /h. This 
property is related to the "double step" of the usual Hall 
conductivity at zero density (electroneutrality) , which is 
a hallmark of the quantum Hall effect in a symmetric 
graphene bilayerpJ 

Another quantity we can easily compute is the total 
magnetization M. of the K + valley. Recall that at fi- 
nite U each state {s 1 , s 2 , fe} carries the orbital magnetic 
moment M s However, when computing the val- 
ley magnetization at given fixed /i, one must account for 
the Berry phase, which effectively modifies the density of 
states. The net result is that, in addition to summing the 
magnetic moment over the occupied states of the origi- 
nal spectrum, there is an ad ditional contribution related 
to the Berry curvature! 34 1 77 1 Namely, M = M M + Mn, 



where 



M 



sr f d k 

9s ^ J J2nji 



M a OU-E a 



(84a) 



a = {s 1; s 2 , fe} is a composite index, 
d 2 k 



Sc ^ 7(2l)2 ("--BJ + ^. (Mb) 



and F + = FQ(F). Thus, for fi > U, in which case the 
occupied states of both conductions bands fill a circle, 
the integration limits are from k = to k = fc F s where 



k^ 2 + U2 - S ^ ] 



(85) 



For E+ < fi < U , the occupied states fill an annulus in 
momentum space. The limits on k are from the inner 
radius k F , to the outer one k F _ . 
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Using the relation 



A4q to the integral over the Berry phase: 



n 



1 d$ 



B 



2itk dk ' 



(86) 



which follows from Eq. ( 79 ) , we reduce the expression for 
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FIG. 13. (Color online) At energy \y\ > Ej, there can be one 
or two branches of k(p), each with its own (a) contribution 
to the anomalous Hall conductivity. As a result the (b) total 
anomalous Hall conductivity has cusps where Fermi surface's 
topology changes. The same also occurs for the (c) Magneti- 
zation density M (fj.) as function of chemical potential fj,. The 
solid curves in each are for U w O.leV. a H for U ~ 1 meV 
is also shown in (b) as dashed trace that it approaches the 
universal value of 2. The negative value for A4(0) is due to 
contribution from the two filled valence bands. When contri- 
bution from all four bands are included, M eventually adds 
to zero. 



M = 



he 



E 

s l i s 2 



(2tt) 2 k 2tt ^ 



(87) 



At this point we recall that the orbital magnetic moment 
M„ given by is related to the difference of the semiclas- 
sical and Berry phases, see Eq. (54). As a result, the 
desired combination M M + M. n is given by the integral 
over the semiclassical phase: 

d 2 k 



M 



9s? 
he 



E 



{2-Kf 



2tt 



Q(ti-E a 



which can be evaluated in closed form. The contribu- 
tion from the two (partially occupied) conduction bands, 
using g s = 2, is 



M(fi,U) 



eJl + _e_ £^ 
ttHc ttHc r y 1 

(0 



x < 



2^? 



El 



El 



'7i 



2UE* 
7i 



(89) 

U < n < E , 
E <^. 



Here the first term, which is linear in [/, is due to the fully 
occupied valence bands, while the additional four possible 
contributions describe the contribution of the conduction 
bands. Interestingly, once the higher energy band s 2 = 
+1 becomes occupied, fi > E^, the total magnetization 
no longer depends on /i, because of partial cancellation 
between the two conduction bands. 

The function M(fJ,,U) at U = 80mcV is plotted in 
Fig. [l3|(c). Similar to the Hall conductivity, it has cusps 
at the energies where the Fermi surface topology changes. 
Specifically, for the first two of them we find 



eU 

TtHc 



M{U) 



eU 



irhc + 4U 2 



(90) 



In Fig. |13[c), the sign of M is negative. However, this 
is unrelated to either paramagnetism or diamagnetism 
because the external magnetic field is assumed to be zero, 
in which case the K~ valley makes an equal and oppo- 
site contribution to the total magnetization of the sys- 
tem. Only the square of M a contributes to the magnetic 
susceptibility: 



Xp = M 2 a v= ^g 2 ii 2 B v. 



(91) 



But \p is only one of the terms (known as the Pauli para- 
magnetism) which dete rmine magnetic susceptibility. As 
shown in previous workj 35 ! 81 " 83 ! the total susceptibility x 
of BLG also contains the Landau diamagnetic term 



XL 



(m e /m e g) h%l> , 



(92) 
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as well as other contributions, which together generate a 
very complicated dependence of \ on (Here v is the 
total electron density of states at the Fermi energy and 
l/ m cff — dE^/dk 2 is the inverse effective mass.) 

Concluding this section, we note that M a , ct h (/j,), and 
M{fJ-) in BLG were previously calculated numerically in 
Ref. [3U Our analytic results for M a agree with that 
work. (For the ease of comparison, a second axis is in- 
cluded in Fig. [8]) On the other hand, there are noticeable 
differences for c H (/i) and M(fi). Regarding c H (^), we 
suspect that the authors of Ref. [34] included the effect of 
impurity scattering in the form of the side-jump, which 
we ignore. The plot of M (^i) presented in Ref. [M] lacks 
the cusps that should be there due to the changes in the 
Fermi surface topology, see our Fig. |l3|c). 



VI. DISCUSSION AND CONCLUSIONS 

In this paper we have presented a quasiclassical Lan- 
dau quantization procedure which includes both the 
Berry phase and the magnetoelectric effects on the band 
structure. This method provides an intuitive picture of 
the Landau level dispersion and several other measurable 
properties of biased BLG. In some cases, we have been 
able to derive analytic expressions for the Landau level 
energies; we also discussed how they may be computed 
numerically. 

Our results are applicable in the analysis of a number of 
experiments which probe transport and thermodynamic 
properties of BLG, including cyclotron resonance, acti- 
vated conductivity, charge compressibility, and magneti- 
zation. Of course, a more realistic calculation of these 
quantities should also include interaction effects. The 
self-consistent mean field approximation for BLG has 
been addressed in several published works, but gener- 
ally such treatments have neglected exchange and corre- 
lation effects, which were considered in Refs. [84"H8T1 and 
shown to give as much as a ~ 30% correction to the mean 
field (Hartree) approximation, similar to the case in two- 
dimensional (2D) electron systems in semiconductors^^ 
Currently, experimental results for the Landau level en- 
ergies from the cyclotron resonance^ and the charge 
compressibility studies^ can be fitted to the theory if 
undetermined variables (U, for example) are treated as 
adjustable parameters. Incorporating all major exper- 
imentally relevant ingredients - Hartree, exchange, and 
disorder contributions — into the same calculation would 
be a more stringent test of the theory. 

Although the Landau level dispersion and therefore 
Landau level crossing points cannot yet be calculated 
with a high degree of accuracy, phenomena that may 
be observed at such points are quite interesting. Indeed, 
crossing of Landau levels has been previously studiecP^2l 
in the context of the quantum Hall effect in conventional 
2D systems. In those systems, the crossings are between 
Landau levels of different subbands or between spin-split 
levels of different Landau levels of the same subband. 



Near the crossing the energy gap vanishes, and so a spike 
in conductance is expected. In the quantum Hall effect 
conditions, this is simultaneously a spike in resistance. In 
experiments, such spikes have been observed to be hys- 
teretic. Sometimes, they were also accompanied by a spa- 
tial anisotropy of the transport. The leading theoretical 
explanation"^ 3 - attributes these phenomena to quantum 
Hall ferromagnetism (QHF). Namely, when two Landau 
level are nearly degenerate and the chemical potential is 
close to the crossing energy, the occupation of the Landau 
levels are modeled as two states pseudo-spin system. De- 
pending on the nature of the crossing, QHF can be of ei- 
ther easy-axis or easy-plane type. In the former case, one 
expects formation of domains whose collective dynamics 
can in principle generate both hysteresis and anisotropy. 
The BLG appears a promising system to study QHF be- 
cause of its high tunability and a rich pattern of level 
crossings we discussed in the paper. 

We are particularly grateful to F. Guinea for valuable 
interactions in the early stages of this work. We also 
thank V. Fal'ko and Q. Niu for discussions. This work is 
supported by the NSF under Grant DMR-0706654 and 
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Appendix A: Low-energy theory of BLG 

In this section we derive the low-energy of BLG by the 
standard method of canonical transformation. Our re- 
sults are in a good agreement with previous workPEHHl 
Some minor discrepancies can be attributed to typo- 
graphic errors therein or differences in notations. 

We begin with the Hamiltonian of Eq. ([3|. The bi- 
layer's electronic structure has four bands. When q = 
K ± , the two central levels lie at E — ±U. For \S q \ <C 1, 
where Sq is given in Eq. Q, we can derive an effective 
2x2 Hamiltonian by writing 



H K+k - H 



V. 



(Al) 



where H° = H(K ± ) contains terms dependent on r y 1 , A', 
and U, and V contains the 7 <S'<j, J^S q , and terms. 
(To lighten notations, the subscript q is dropped in the 
following.) 

The unperturbed Hamiltonian H° has levels E 4 4 = 

T\/7i + U 2 and E% 3 = A'^fU. The eigenf unctions \ip,) 
are the column vectors of the matrix 



f = 



10 

cos(0/2) sin(0/2) 

-sin(0/2) cos(0/2) 

10 



(A2) 



where tan$ = Ji/U. Eliminating the high energy sub- 
space spanned by V*! 4) by unitary transformation 



H = e lQ He~ lQ . 



(A3) 
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we obtain the effective 2x2 Hamiltonian 

TJ 77*0 C I T/ 



+ X/ ( RO _ 



1 



(A4) 



n a 



up to terms of order V 2 . Here n,n' G {2,3} are labels 
for the low energy subspace, while a € {1,4} labels the 
high energy subspace, and = (V^ |V| ipj). The matrix 
elements of Q are given by 



Qua = 1 — n Sn + * 2^ 



' 51 i c?o 



(A5) 



/(£°-£°)(i?° -£0) 



(V 2 ) , 



with Q a „ = (Q on ) f ; 

Proceeding in this manner, we obtain the 2x2 block 
for the inner bands, 



H = 



e - U + uj 



r e Q + U-uj 
where, to lowest order in U and A', 



(A6) 



I 7 ^ + (7 ° + J ] A ' ) {s, } + u -4 [S, gt] , 



7i 



7? 



ll 



{S, St} 



(A7) 



7i 



27f 



V 7i 7i / 



(A8) 
(A9) 



Anticipating the introduction of an external magnetic 
field, we have allowed for the possibility that S and S> 
do not commute. Recognizing that 74/70 = 0.05 <C 1, it 
is permissible to drop the terms of order 7! and 74 A', in 
which case 



£n = 



7o_A 
27? 

ll 



{S,S^} 

{s, st}- 



^71 
tIa 

27i 2 



[5, ^] 

[s, st] 



^ = 7 3 St-is 2 , 
7i 



(A10) 
(All) 
(A12) 



leading to Eq. ( 33 1 . Our results agree with those of Ref. 12! 



if A' and 74 are set to zero. 

For B = 0, in the vicinity of the K ± points, the four 
bands disperse as shown in Fig. [I] The two central bands, 
which comprise the low energy sector, are separated by 



2U at fc = 0. Their dispersion is described by the effec- 
tive Hamiltonian of Eq. ( A6 1 . One finds that for fc = k x 
the central bands have a characteristic double hump (or 
Mexican hat) shape provided U > 27-^3 « 80meV. 

It is convenient to write H = e Q + 23(fc) cr, where B z = 
lu, B x — iB y = £, and cr is the vector of Pauli matrices. 
When the actual magnetic field B vanishes, 



_ a 4 

lo 



if 

ll 



(A13) 
(A14) 



where e k = Hvok, the origin in fc-space is taken as one of 
the K points, and ip = tan -1 (fc / k x ) is the correspond- 
ing polar angle. The eigenvalues are 



E. 



Sl \B(k) 



(A15) 



Let us now discuss the Berry phase. Semiclassically, in 
the presence of a weak magnetic field, the wavevector fc 
evolves in time according to Eq. ( 59 1 : 



fc = LO r z x fc 



2tt 



sgn(v 



If we can neglect 73, then the trajectory the pseudospin 
traces on the Bloch sphere winds twice for every cycle 
of fc, owing to the e 2lip factor in B x — iB y = £. There- 
fore, the accumulated Berry phase is equal to 2 X 4 = 1 
times the solid angle traced by vector 23 (fc). Actually, 
the Berry phase is defined modulo 2ir. To be consistent 
with the earlier choice of the overall phase factor of the 



basis state (49), we need to subtract 2ir from the solid 



angle. The result is 



$b = 2tt( 1 



1251 



2tt 



2ttU 



E 



M 
7? 



(A16) 



However, it differs from our earlier Eq. (67) for the Berry 



phase in BLG. The discrepancy arises due to the canoni- 
cal transformation by which we obtain the wavefunction 
in the new basis: \tp') = er tC! 



(5$ = $ 



2tt- 



el B r 



B 



7i 



\B\ 



—2iri (t/)' |e 



iQ 



-%Q\ 



(A17) 



The combined phase $ B = $^ + <5$ is in agreement with 



the four-band expression Eq. (67), to within the accuracy 
of this calculation. 

Finally, we can go beyond the semiclassical approxi- 
mation, obtaining the effective Hamiltonian in Eq. (36). 



If 73 is neglected, each pseudospin component Landau 
level is connected to a unique mate, and the Hamilto- 
nian breaks up into a direct sum of 2 x 2 blocks, given 
by Eq. (38) (n = — 1 and n = are special cases where 
H reduces to a scalar). 
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FIG. 14. (Color online) Sketch of the structure of the mag- 
netic bilayer Hamiltonian, showing nonzero matrix elements 
as links. Each link between orbitals in column n and column 
n + 1 is multiplied by a factor \Jn + 1 (a>o/7o)- The diagonal 
entries in the Hamiltonian for each orbital are given at the 
upper left. When 73 = 0, the Hamiltonian breaks up into a 
direct sum of 4 x 4 blocks. 



Appendix B: Matrix representation of the 
Hamiltonian in a finite magnetic field 

In the presence of a magnetic field, the full Hamiltonian 
H + in the K + valley is given by Eq. (16), 



n- 



I —U —uj a i] 4 uj Q a (] 3 w aA 
-LJ al —U + A' 7j ?74W a 
f -v. U + A' -uj n a 



?7 4 w o' 
V % w o a 



7i 
774 w a 



w . 
E7 



/ 



where a and are Landau level lowering and raising op- 
erators, respectively. In the occupation number basis |n), 
the matrix elements of H + can be understood pictorially, 
by referring to Fig. |14| Writing the general wavefunction 



I*) 



E 

n=0 



( u n \ n ) , V n \n) , u n \n) , v n \n)) T , (Bl) 



the links in Fig. [14] indicate matrix elements between the 
various components {u n , v n , u n , v n }. 

One finds that W. = W. a © 'H h © T-L c can be written 
as a direct sum of three terms. In evaluating the spec- 
trum numerically, we truncate % ai b,c at a high Landau 
level index, as shown in the figure. Typically we chose a 
maximum index of n max 300, checking that spectrum 
did not vary significantly as the upper index cutoff was 
further increased. This feature is most evident at high 



FIG. 15. (Color online) Landau level energies vs. interlayer 
bias U for a field value B — 20 T. Solid lines correspond to 
the K + valley and broken lines to the K~ valley. The color 
and line type are as in Fig. [6] The shaded area indicates the 
energy gap at B = 0. 



fields, such as in Fig. [T5j where we have taken B = 20 T. 
The spectrum of "H a is shown in black, that of "Hb in red, 
and that of TL C in blue. Solid lines correspond to the K + 
valley and broken lines to the K~ valley. One sees in the 
figure that curves of the same color and line type cannot 
cross at an accidental degeneracy. 

Frequently in this paper we have ignored the SWMc 
parameter 73, setting it to zero, In this approximation, 



as can be seen from Fig. 14 the occupation number space 
Hamiltonian further resolves itself into a direct sum of 
4x4 blocks, given by the expression in Eq. (27), which 
connect {u n _ 1 , v n , u n , v n+1 } for each n. (There is also 
a remaining lxl and 3x3 block associated with the 
indices n = and n = 1.) 

In addition to eigenvalues, we also calculated the eigen- 
functions, one of which is shown in Fig. |11| To do so we 
chose the symmetric gauge, where the Bloch wavefunc- 
tions of |m) oscillator states are given by 



nm+l 



lk\n) 



V2 



-fc 2 £|/4 



(B2) 



These basis states were weighted with the coefficients ob- 
tained from diagonalizing the Hamiltonian matrix and 
then summed over all components (both n and the sub- 
lattice index). 

From these calculations we concluded that the effect of 
7 3 diminishes as B increases, as was previously observed 
in Ref. [2] The semiclassical argument that explains this 



behavior was given in Sec. IV B Here we mention another 
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reasoning P which is based on the usual perturbation the- 
ory. 

The leading-order correction to the energies due to 7 3 
is approximately u>q/ AE, where AE ss fyv g \/(ktfg) is the 
Landau level spacing. Therefore, the relative magnitude 
of this energy shift is small provided 



73 « 



7o 



(B3) 



737i/7oi which is roughly consistent with the threshold 
B ~ 1 T where the effect of 7 3 is observed to become 
insignificant in the numerical calculations. On the other 
hand, at finite U and near the bottom of the Mexican 
hat, where v g — 0, the expression on the right-hand side 
of Eq. (B3| diverges. This implies that the effect of 7 3 



is larger and persists to higher B. This is also consistent 
with the numerics, see Sec. |IVB| 



At U — and e k <C "f 1 this inequality gives 2 wq 3> 
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